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Abstract Prior to the detection of its outermost Uranus-mass object, it had been 
suggested that GJ 876 could host an Earth-sized planet in a 15-day orbit. Observation, 
however, did not support this idea, but instead revealed evidence for the existence of 
a larger body in a ~125-day orbit, near a three-body resonance with the two giant 
planets of this system. In this paper, we present a detailed analysis of the dynamics 
of the four-planet system of GJ 876, and examine the possibility of the existence of 
other planetary objects interior to its outermost body. We have developed a numerical 
scheme that enables us to search the orbital parameter-space very effectively and, in 
a short time, identify regions where an object may be stable. We present details of 
this integration method and discuss its application to the GJ 876 four-planet system. 
The results of our initial analysis suggested possible stable orbits at regions exterior 
to the orbit of the outermost planet and also indicated that an island of stability may 
exist in and around the 15-day orbit. However, examining the long-term stability of 
an object in that region by direct integration revealed that the 15-day orbit becomes 
unstable and that the system of GJ 876 is most likely dynamically full. We present the 
results of our study and discuss their implications for the formation and final orbital 
architecture of this system. 

Keywords Stability ■ Resonance • Hamiltonian Systems • Numerical Methods • 
Planetary Systems 



1 Introduction 



Since the announcement of its two resonant planets bv lMarcv et"al~l l|200lh , the plan 



etary system of GJ 876 has had a special place in exoplanetary science. As the first 
system detected with two planets in a mean- motion resonance (MMR), [planets GJ 
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876 b and c, with minimum masses of 2.27 and 0.71 Jupiter -masses and orbital periods 
of 61.11 and 30.08 days, respectively ( Rivera et al. I . lioiol ). are in a 2:1 mean-motion 
resonance] , G J 876 has been the subject of extensive research and has had major con- 
tributions to the development of the m odels of planet migration and resonance capture 
( Lee fc Peale]l2002l . iKlev et al. If2005h . Because of the short period gravitational inter- 
actions between its two giant planets, the dynamics of this system and the possibility 
of its hosting additional planetary bodies ha ve been subjects of i ntense studies. Soon 
after t he announcement of its 30-day planet ( Marcv et al. l . fioOll 'l. lKinoshita fc Nakai I 
(|200lh studied the stability of this system and showed that its 2:1 MMR presents the 
most stable orbital configuration for different values of its planetary masses and orbital 
elements at different epochs. .Rivera fc Lissa uer (200ll) andlLaughl in fc Cham bers I 
( 200 ih also studied the dynamics of this system. These authors argued that because of 
the relatively small separation between the two giant planets and the possibility of their 
close approaches, when fitting the system's radial velocities, the interactions between 
these two planets ha ve to be taken into accou nt. By presenting new fits to the radial 
velocities of GJ 876, 1 Rivera fc Lissauer ] (|200ll ) showed that the system will be stable 
for long times when the planet s' orbits are co-planar and ha ve low inclinations with 



respect to the plane of the sky. iLaueh in fc Chambers J200lh Umited this inclination 



to a range of 30° to 53°. The results bv lRiverafcTissauer I a rid Laughlin fc Chambers I 
were later confirmed by the analytical mo dels of Ji_et_aL (l2002|') an d in a thorough 



study of the global dynamics of GJ 876 bv lCozdziewski et al. l|2002l '). 

Being a dynamically interesti ng system, GJ 8 76 has also been the target of ob- 
servation for many years. In 2005, iLaughlin et al. I re-analyzed the new radial velocity 
data of this star and showed that when planet-planet interaction, stellar jitter, and 
instrumental uncertainties are taken into account, a 2:1 resonant co-planar orbit with 
an inclination less than 20° would present the most viable and stable planetary con- 
figuration for the two giant planets of this system. These authors also showed that in 
addition to being in a 2:1 MMR, the two planets are locked in a secular resonance where 
they librate around apsidal alignment with an amplitude of 34°, and their joint line of 
apsides precesses at a rate of 41° ev ery year. The existence o f this s e cular resonance had 
also b een p resented in the work s of lLaughlin fc Chambers"! ( 200lh . [Gozdziewski et al~l 
(I2OO2I ). and I Lee fc Pealel ( 2002h where these authors analyzed the dynamics of GJ 876 
and presented a model for the migration and resonance capture of its two giant planets. 

The con tinuous observ ation of GJ 876 resulted in even more fundamental discover- 
ies. In 2005. [Rivera et al~l announced that a 6.83 Earth-masses planet exists in a 1.94- 
day orbit around this star. This discovery that marked the detection of the first super- 
Earth planet, prompted many astrono mers to examine t he possibility of the existence of 



other Earth-like bodies i n this system dji fc Liu ll2006l. 20071 : iRivera fc Haghighipour I . 

In that direction, iRivera fc Haghighipour I (|2007l ) studied the dynamics of ficti- 
tious planets in the system of GJ 876 and showed that in addition to a small region 
interior to the orbit of planet c where the super-Earth planet of this system has a stable 
orbit, a region of stability exists beyond the orbit of the outer giant planet (i.e., planet 
b) corresponding to an exte rior 2:1 resonance with this object. In 2010, the prediction 
bv lRive ra & Haghi ghipour I materialized and a new Uranus- mass planet was discovered 
in a ~ 125-day orbit around GJ 876, making this system the first planetary system 
with three planets in a Laplace resonance. 

Prior t o the detection of its fou rth planet, the three-plan e t system of GJ 87 6 was 
studied bv lBean fc Seifahrt I (|2009l ') and lCorreia et al. I j 20 id ). I Bean fc Seifahrti stud- 
ied the architecture of the system and showed that the assumption of co-planarity is 
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valid, and the mutual inclination of the two giant pla nets axe within 1° to 7°. In a 
very thorough dynamical analysis, ICorreia et al. 1 (|2O10h used the combined data from 
HARPS and Keck, and while confirming the existence of planet d in a 1.94-day orbit, 
showed that the orbits of the two giant planets of the system are co-planar with incli- 
nations of approximately = 48.9° and ic = 48.1°. These authors also suggested that 
an Earth-sized planet in a 15-day orbit at ~ 0.08 AU could have a stable orbit for a 
long time. This finding was very interesting since the existence of such an additional 
planet would help to explain the anom alous high ecce ntricity of planet d. We would 
like to mention that in a recent article, iBaluev | (|201lh re-analyzed the HARPS data 
and suggested that the non-zero eccentricity of planet d may be due to the lack of 
proper interpretation for red-noise in t he data. 

With their new four-planet system, iRivera et_aL_ (l2O10h examine d the possibility 
of a stable solution for the planet proposed bv ICorreiaet al. I ( 2010l ). However, they 
were unable to find a signal corresponding to the 15-day planet in their radial velocity 
observations. Placing test particles in the region around 0.08 AU and integrating the 
entire system numerically, these authors found that only a small fraction of particles 
remained stable for the 10 Myr duration of the integration. 

Given the resonant state of the three planets of the system, and that these planets 
have most probably captured each other in resonance while migrating inwards from 
outer regions, and also given the fact that GJ 876 is an M star and planet formation 
around M stars favors the formation of low-mass objects, it would be natural to examine 
whether GJ 876 can host additional (low-mass) planets, in particular an Earth-like 
object in a 15-day orbit. This paper addresses this and several other questions regarding 
the dynamics of the planetary system of GJ 876. 

In the rest of this paper, we present our methodology and the results of our extensive 
numerical analysis of the dynamics of the system. We are particularly int erested in the 
dynam ical surrounding of the fitted orbit of planet e as given in Table 3 of lRivera et al. I 
(2010), and the pos sibility of addit i onal p lanets in the system including the Earth-sized 
planet proposed bv lCorreia et al. ] (:2010). We discuss our numerical method in section 
[21 and present the results in section [3] Section |3] summarizes our analysis where we 
also discuss their implications for the formation of this planetary system. 



2 Numerical tools - Stability maps from variational equations 

As mentioned above, the goal of our study is to examine the possibility of the exis- 
tence of additional bodies in the four-planet system of GJ 876. Our approach is to use 
stability analysis to identify regions where an object can have a long-term stable orbit. 
Since there is no analytical solution to a general A^-body problem with TV > 2, the 
orbital evolution of the planetary system of GJ 876 has to be computed numerically. 
As explained below, we will use a symplectic integrator for this purpose. Due to their 
reliability and stability properties, especially for long-term integrations, these integra- 
tors have fou nd a special place in celestial mechanics and planetary dynamics. We refer 



the reader to lHairer et al. 1 (l2002h for a gener al overview of symplecti c integrators and 
to, for instance. iGladmanetaT' ( 199 ih . and lEggl fc Dvorak! (|2010l ') for applications 



of these integrators to celestial mechanics. 

To study the dynamical stability of an object, one has to identify (and exclude) 
conditions that result in chaotic motion of the body. Several methods exist for this pur- 
pose that discriminate between regular and chaotic motions. One can either analyze 
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the o rbit of a particular object [using e.g., the frequency analysis technique (|Laskar I . 
I1993I I]. or study the evolution of deviation vectors (S) from a given orbit. Examples of 
the l atter inc l ude th e maximal Lyapunov Characteristic Exponent (7) [for a review, see 
e. g. ISkokos I 1 201Ci )]. the Fast Lyapunov Indicator te chnique ( FLI) ( Froeschle et al.~l . 
Il997l '). the Sm aller Ahgnmeiit Inde x (SALI) method ( Skokos I . boOlh and its general- 
ization GALI (jSkokos et al. 1 . 120071') . and the MEGNO (M ean Exponential Growth of 
Nearby Orbits) chaos indicator ( Cincotta &i SimdT l20od ) . 

In this paper we conc entrate on MEGNO as a fast and reliable method to study the 
dynamics of single orbits ( Mafflone et al. Ll201lh . For this chaos indicator, the variation 
of a quantity Y and its mean value Y during time t are given by 



■ s ds 



and 



m = - 



Y{s)ds. 



(1) 



10 " ' " '-Jo 

The quantity S in equations ([T} denotes the derivative of the deviation vector with 
respect to the independent variable s. To determine the stability of an orbit, one has 
to analyze the time evolution of Y , which is connected to the maximal Lyapunov Char- 
acteristic Exponent of the orbit. For stable, periodic orbits, Y approaches asymptot- 
ically, while for quasi-periodic ones it will tend to 2. For chaotic initial conditions, Y 
will grow with time as 'yt/2. The Lyapunov time, defined as Tj^ — 7^^, can therefore 
be obtained from Y using 



Tl 



t 

W' 



(2) 



In this equation, t is the total integration time. 

The above-mentioned variational equations have to be integrated together with the 
equations of motion of the system. To increase the speed of these computations, we 
used a sympiec tic integrator, and integra ted the entire set of equations of the system 
simultaneously. Tskokos fc Gerlachl (2010|) discuss different methods of employing these 
integrators for advancing the variational equations of a Hamiltonian system. The most 
ef ficient of these methods, the s o called Tangent Map (T M) technique, was applied 
bv lSkokos fc Gerlach I ( 20ld l and lCerlach fc Skokos I ( 2011) mainly to low - dimen sional 
Hamiltonian systems with 2 and 3 degrees of freedom. Gerlach et al. I ( 20111 ) have 
shown that this method can also be applied to multi-dimensional Hamiltonian systems 
and is very efficient and superior to other commonly used numerical schemes, both in 
accuracy and speed. We not e that the TM m e thod was first discussed in the context 
of celestial mechanics by .Mikkola fc Innanen ] (|l999l ). 

In this paper, we explore the dynamical stability of hypothetical planets in the sys- 
tem of GJ 876 for a wide range of their orbital elements. To make thi s task computation- 
ally fe asible, we use the SABAn/SBABn integrators developed bv lLaskar fc Robutel I 



(|2001h . which proved to be efficient and reliable. These sympiectic splitting integra- 
tors have been designed specifically for integrating Hamiltonian systems of the form 
H = A -\- eB, where A and B are separately integrable, and e 1, as in the case 
of hierarchical A'^-body systems. For more details on this technique and on different 
me thods of splitting the Hamiltonian into two integrable parts, we refe r the reader to 
e.g. lDuncan et al. I j 19981) . IChamberTI j 19991) . ICozdziewski et al. I (|2008l ) and references 
therein. 

As mentioned above, we used the TM method to compute the variational equations. 
The formulas fo r advancing variational equ atio ns using a sympiect i c inte grator can 



ihe tormulas to r advancmg variational equ atio ns usmg a sympiect i c mte grator can 
be found in e.g. iMikkola fc Innanen I l| 19991 ) and ICozdziewski et al. I 1 20081 ). We used 



5 



these formalisms, specifically the latter, where also a time-discrete approximation of 
equations ([T]) is given. Let us finally note that symplectic methods cannot be used with 
a trivial automated step-size control. Therefore, they are usually implemented with a 
fixed integration step, which is denoted by r in this paper. 



3 Stability analysis 

In this section, we present the results of our stability analysis for the GJ 876 p lanetary 
system. We used the SABA4 symplectic integrator ( Laskar fc Robutel I . I2OO1I ) to inte- 
grate the equations of motion and the variational equations of the system. The latter 
are computed only for the test particles, which we define in this study to be always 
massless and used to determine the stability in the specific regions. Note that all orbital 
elements are given with respect to the central body and integrations were carried out 
incorporating only Newtonian gravity. No relativistic and tidal effects were included. 



3.1 Global stability in the inner three-planet system 
of GJ 876 

O ne of the exciting results of the dynamical analysis of GJ 876 bv lCorreia et al.~l 
( 2OIOI ) is the suggestion that a stable region exists for an Earth-sized planet at ~ 0.083 



AU. This island of stability that corresponds to an orbital period of 15 days, would 
allow a terrestrial planet to be in a three-body 1:2:4 resonance with the two giant 
planets of the system. The interesting fact is that if this 15-day planet exists, it may 
be possible to use the interaction between this object and the innermost planet of the 
system (the 6.5 Earth-masses super-Earth that orbits the central star in a 1.94-day 
orbit) as a way to account for the relatively high ec centricity of the latter body. 

Despite the predictions of ICorreia et al. 1 ( 201G| ). radial velocity obser vations have 



not be en able to detect a planet in the 15-day orbit. Instead, observations bv lRivera et al. 
( 2010l ) revealed that an additional planet does exist in the system, but it is Uranus- 



mass and in an orbit with a period of approximately 125 days. This body is cu rrently 
the outermost planet of GJ 876 system, and as shown bv iRivera et al. 

] l|2010l l. forms 

a three-body Laplace resonance with the two inner giant planets. 

The purpose of our study is to determine how the discovery of this new planet affects 
the stability of the 15-day orbit. In other words, we would like to examine whether an 
Earth-sized planet can maintain a stable orbit a t 0.08 3 AU in this new four-planet 
configuration? A few simulations bv lRivera et al. 1 1 2010l ) suggested that stable motion 



may be possible around ~ 0.083 AU. However, those simulations had been carried out 
only for short integration times. Our goal is to study also the long-term stability of 
orbits in this region. 

As we mentioned in section [J] our approach is numerical. In order to verify that 
our numerical method produces reliable results, we first used our code and integrated 
planets b, c and d together with a large battery of test particles in the region interior 
to the orbit of the 30-day giant planet. We carried out two se ts of simulations: one 
with the orbital elements for the planets as gi ven in Table 2 of Correia et al. I (|20ld ) 



(hereafter model I), and one with the data from lRivera et a~ ] l|20ld . Table 2) (hereafter 



model II). The integrations were done for t = 500 years using a time-step of r = 0.05. 
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Fig. 1 Stability graphs of test particles in the three-planet system of GJ 876. The color-coding 
represents the mean value of a particle's MEGNO. Red corresponds to chaotic motions and 
black/blue denote stable regions. The integrations were carried out for 500 years using a step- 
size of T = 0.05 days. The initial orbital elements of the planets are given in Tabled Graph (a) 
shows the results for model I and graph (b) depicts those of model II. A grid of 600 x 100 test 
particles was used, for which the initial inclinations were set to j = 50° in the simulations of 
graph (a) and they were considered to be co-planar with the known planets of the system, i. e. 
i = 59°, in the simulations of graph (b). All other orbital angles were initially set to zero. The 
white lines in both graphs mark the boundaries of orbit-crossing with the 30-day and 1.94-day 
planets of the system. 



Figure [T] shows the resuHs. In this figure, the value of the mean MEGNO associated 
with the orbit of each test particle is shown for different values of its semimajor axis 
and eccentricity. As explained before, orbits with low mean MEGNOs are considered to 
be stable whereas those with high values (e.g., red color) are irregular and/or chaotic. 
The left panel of Fig.[T]shows the results for the three-planet system of model I and the 
right panel shows similar calculations for model II. The initial orbital inclinations of 
test particles in the left panel were set to 50° whereas in the right panel test particles 
were considered in the same plane as the planets of the system, i. e. i = 59° . All other 
angular elements of these objects (i.e., longitude of periastron, argument of ascending 
node, and mean-anomaly) were considered to be zero. 

As shown by the two panels of Fig. [TJ an island of stability exists in the region close 
to 0.08 AU. While in both simulations, the width of this region is almost identical, it 
seems that in the simulations of Fig.[TJa), where the orbital elements of the planets were 



Table 1 Orbital parameters for the different models used in sections 13. II and 13.21 

model I model II 

Correia et al. ('2010. Table 2) Rivera et al. (2010. 



M* [Msun] 




0.334 






0.32 




planets 


b 


c 


d 


b 


c 


d 


a [AU] 


0.211 


0.132 


0.021 


0.208 


0.130 


0.021 


e 


0.029 


0.266 


0.139 


0.029 


0.255 


0.257 




48.93 


48.07 


50 


59 


59 


59 







-2.32 














a;[°] 


275.52 


275.26 


170.60 


50.70 


48.67 


229.0 


A n 


35.61 


158.62 


29.94 


16.10 


343.67 


227.0 


M [Mj,p] 


2.64 


0.83 


0.0198 


2.276 


0.714 


0.022 



7 



taken from lCorreia et al. I (|2O10l '). the stability does not include orbits with eccentricities 
close to zero. In Fig. [TJb), however, we find that the most stable test particles in the 
region of 0.08 AU are the ones in circular orbits. The apparent instability for zero- 
eccentricity orbits in Fig.[TJa) could be attributed to the differences between the orbital 
elements used in the two models. For instance, while in model I a full three dimensional 
solution for the planets of the system is used, in model II, it is assumed that all planets 
are on the same plane. These differences could cause slight variations in their regions 
of stability. 

We recall that the purpose of carrying out the above-mentioned simulations and 
generating Fig. [2a) was to test the reliabili ty of our num erical method. A comparison 
between this figure and Fig. 10 of .Correia et al. indicates that our integrator 

is in fact reliable as it has been able to correctly re-produce the main features (i.e., the 
width and location of the island of stability at 0.08 AU) of the figure published by these 
authors. The slight differences between the two results, such as the stochastic pattern 
at a « 0.05 AU, are primarily due to the intrinsic chaotic characterist ic of the problem 
and the differences in the used integrating schemes. As shown bv ICerlach I (|2008l ). 
already existing small differences, as those resulting from the use of different chaos 
indication methods, or different numerical schemes, etc., could lead to very different 
stability results in areas where the phase-space is highly structured. 



3.2 Orbital stability in the outer three-planet system of GJ 876 
As indicated bv iRivera et al. 

] (|2O10l '). the orbit of the recently detected Uranus-mass 
planet of the system is stable for at least 400 Myr. Test particle simulations led these 
authors conclude that in order for this planet to maintain orbital stability, its orbit 
has to be in a Laplacian resonance with the two giant planets of the system where its 
period will be in a 4:2:1 ratio with the orbital periods of these bodies. 

To verify the resonant state of the new planet (hereafter, planet e), and to quantify 
the size of its stable region, we integrated together with the massive planets of model 
II, the orbits of a large number of test particles in the region exterior to the orbit of 
planet b (orbital period of ~61 days) for different values of their semimajor axes (a) 
and mean longitudes (A). We considered the test particles to be initially in circular 
orbits, and co-planar with the 3 planets (i. e. i — 59°). We set all other orbital angular 
elements equal to zero. We integrate the system for 10^ orbital periods of the proposed 
planet e (« 3400 years) using a step-size of r = 0.1 days. 

Figure [2] shows the results. The (a. A) stability map in this figure clearly shows the 
boundary between regular and chaotic orbits. Using equation 0, one can estimate 
that the stable test particles in this figure have Lyapunov times of > 340 years. 
An interesting feature of Fig. [2] is the small island of stability around a = 0.33 AU 
where planet e exists. The value of A for this region is ~ 220° which places it in a 
Laplace resonance with planets b and c. Given the small size of this area, it is clear 
that this resonance protects planet e against close encounters with the other planets 
of the system. 

Although in their analysis, I Correia et al. ] (|2010l ) did not discuss the possibility of 



an island of stability at 0.33 AU, one can see from their Fig. 10 that the simulations 
by these authors also show a small region of stability around that area. The values 
of orbital eccentricity corresponding to this stable region is clos e to zero which agree s 
with the orbital eccentricity of planet e (0.055) as reported bv I Rivera et al. I (|2010l ). 
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Fig. 2 Graph of the orbital stability of test particles exterior to the orbit of planet b. Simula- 
tions were carried out for different values of the particles' semimajor axes and mean longitudes. 
The color-coding corresponds to the particles' mean MEGNO as described in Fig. [T] Integra- 
tions were done for ~ 3400 years using a step-size of r = 0.1 days and a grid of 400 x 150 test 
particles. The initial orbital inclinations of test particles were set to i = 59°, i.e., co-planar with 
the orbits of the three planets. However, all other orbital angular elements were initially set to 
zero. Results show a small stable region at A ~ 220° in a Laplace resonance with planets b and 
c, and in the nominal position of planet e (white cross) as given in Table 3 of iRivera et al. I 
ll201Cl'l . 



This agreement between our results and the results of lCorreia et al. I ( 201Cll l could serve 
as another confirmation of the validity of our integration method and the reliability of 
our results. 



3.3 Could an additional small planet exist in the GJ 876 four-planet system? 

As shown by our stability analysis of test particles in the three-planet system of GJ 
876 (Fig. [T]), an island of stability se ems to exist in and around 0.083 AU, even when 
we use the orbital elements given bv lRivera et al. ] (|2010h . A planet in this region will 



be in or near a 2:1 MMR with planet c, and as a result, will be protected against 
close encounters with other bodies. Given the orbital distance of planet e from this 
stable region, and the fact that the orbit of this object is in a mean-motion resonance 
with planets b and c, it would seem natural to assume that the addition of planet e 
to the system would not alter the stability of a planet at 0.083 AU. To examine this 
assumption, we carried out similar simulations as in section 13.11 but with plan et e 
includ ed. We used the planets' orbital elements as given in Table 3 of Rivera ct al. I 
( 2010l ) (hereafter called model III, see Table and calculated the mean value of 



MEGNO for a suite of test particles in the region between 0.02 AU and 0.1 AU. 
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Table 2 Orbital parameters for the model used in section [331 



model III 

source ^iver^^t_al^ ^2010, Table 3) 

[Msun] 0.32 



planets 


b 


c 


d 


e 


a [AU] 


0.210 


0.130 


0.021 


0.328 


e 


0.037 


0.256 


0.207 


0.045 


i [°] 59 


59 


59 


59 




n [°] 













uj [°] 43.27 


48.74 


234.07 


251.36 




A [°] 15.88 


343.33 


229.38 


214.06 




M [Mjup] 2.276 


0.714 


0.022 


0.046 





Figure [3] shows the results. ^From this figure, it can be seen that the stable island 
at 0.083 AU still exists. However, its size on the eccentricity and semimajor axes has 
become smaller. This figure also shows that the region corresponding to the most stable 
orbits (the darkest area) is now exclusively for test particles in circular motion. Figure 
[3Ib) shows that the stability also crucially depends on the initial value of the mean 
longitude, which should be zero in order to place the body inside the 2:1 MMR with 
planet c. 

Although our calculations of MEGNO suggest stability at 0.083 AU, it is important 
to emphasize that this result has been obtained only from 500 years of simulations. This 
time span is sufficient to identify unstable regions at reasonable CPU costs, especially 
when the grids in semimajor axis, eccentricity, and other orbital elements are very 
dense. However, it is not long enough to allow us to make conclusions about the long- 
term stability of the identified stable regions. For the latter, integrations have to be 




Fig. 3 Graphs of the orbital stability of test particles in the four-planet system of GJ 876. 
Simulations were carried out for different values of the semimajor axis, eccentricity, and mean 
longitude of test particles on a grid of 600 x 100. The color-coding represent the mean value 
of MEGNO as described in Fig. [T] Integrations were carried out for 500 years using a step-size 
of T = 0.05 days. The initial orbital elements of the planets are given in Table [2] All test 
particles were placed in co-planar orbits with planets, i. e. i = 59°. In graph (a), the initial 
values of other angular variables of test particles were also set to zero. However, in graph (b), 
the mean anomalies of the test particles were kept at non-zero values, but instead, their initial 
eccentricities were set to zero. In (a) the white lines mark the boundary, above which collisions 
with planets b and d (crosses) are possible. 
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carried out for longer times. In that respect, Fig. |3]can be used to make the search 
for long-term stability very efficient by allowing to focus only on the areas of the 
orbital parameter-space where an object may have a chance to be stable. For the four- 
planet system of GJ 876, this area is the 15-day orbit. We, therefore, considered a test 
particle with orbital elements corresponding to the most stable zone around 0.083 AU 
(the darkest spot in Fig.fS]), and integrated its orbit along with all other planets of the 
system, for longer times. The results are shown in Fig. 3] by red color. As can be seen 
in this figure, the orbit of the test particle was stable for only 9 x 10^ years. During 
this time, the object was trapped in a 2:1 mean-motion resonance with planet c and 
its corresponding resonant angle, 9 = —A -I- 2Ac — oJc, was librating around 310°. The 
difference between the longitude of the periastron of the particle (u) and that of planet 
c (a;c) at that state was Auj — uj ~ Uc ^ 250°. 

Figure |4] also shows that during the course of integration, strong perturbations 
from planet c caused the eccentricity of the test particle to rapidly increase. At this 
state, the variations in Auj (and 9) developed large amplitudes until shortly before the 
ejection of the particle from the system when Alj and 9 began to circulate. The latter 
indicates that the orbit of the particle was no longer stabilized by the 2:1 MMR. 

As mentioned before, the main purpose for carrying out these simulations was to 
determine the extent to which the stability of an object in the 15-day orbit would be 
affected by the newly discovered planet of the system. In that respect, and in order 
to be able to make a comparison with the stability analysis of the system prior to the 
detection of planet e, it is necessary to carry out similar A^'-body integrations for the 




Fig. 4 Results of the long-term integration of the orbit of a test particle from the most 
stable part of the island of stability in Fig. [S] The semimajor axis of the particle is initially 
a = 0.081420 AU. Its other orbital elements were initially set to zero. The initial orbital 
elements of the four known planets of the system were taken from model III (Table |2]l. The 
integration was carried out using a step-size of t = 0.05 days. The graph on the left shows 
the orbits of the planets and that of the test particle (shown in black) in the x — y plane. The 
graph on the right shows the time evolution of semimajor axis a, eccentricity e, difference in 
periastron longitudes Au} = ui — uic, and resonance angle 9 = —X + 2\c — IjJc of the test particle. 
Shown in black color in the background is the orbital evolution for a test particle with the 
same initial conditions, but integrated in a system with only 3 planets as in the right panel of 

Fig.m 
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15-day orbit test particle with only planets b, c, and d included. We used the orbital 
elements of these planets as in model II and integrated the orbit of the same test 
particle as in Fig. [31 for 1 Myr. These results are presented as black dots in the right 
panel of Fig. [3] We found that for this configuration, the orbit of the test particle was 
stable for the entire duration of the integration. The latter indicates that the addition 
of planet e had a profound effect on the stability of this object. Given that the orbit 
of our test particle was the least chaotic one in the area around 0.083 AU, this result 
further suggests that the region between planets b, c, and e in the system of GJ 876 is 
naturally unstable. In other words, the system of GJ 876 seems to be dynamically full. 

3.4 Stability outside planet e 

We also examined the stability of an object exterior to the orbit of planet e. We used 
the orbital elements of the planets from Table [2] (model III) and carried out similar 
simulations for test particles in the region of 0.2 < a < 2.0 AU. The results are given in 
Fig. [5] As shown there, a large stable region exists outside the orbit of planet e where 
the system can host additional planets. The inner edge of this region is determined 
mainly by interaction with planet e. The fine structures in this region, such as the 
gaps and islands of stability, are due to various mean-motion resonances. For instance, 
the large gap at 0.53 AU corresponds to a 4:1 mean motion resonance with planet b. 
Figure[5]also shows that as the semimajor axis of an object increases, stability extends 
to orbits with larger eccentricities. However, the upper limit of the eccentricity is set 
by the close approach of the object to planet e. 

An interesting feature of Fig. [5] is the small islands of stability that appear at high 
eccentricities outside the gaps of instability that are due to mean-motion resonances. 
Whether these islands of stability are long-term stable is a topic that requires A''-body 
integrations of an actual object in those locations. Such a study is, however, beyond 
the scope of this paper. 

4 Summary and Conclusions 

We presented the results of a detailed study of the stability of the planetary system 
of GJ 876. Our goal was to determine whether an Earth-sized planet could maintain a 
stable orbit in a four-body resonance configuration with the currently known planets of 
the system. We computed the value of the chaos indicator MEGNO for a large number 
of massless test particles to map the system's orbital parameter-space, and identified 
regions where the orbit of an additional object could be stable. Results suggested an 
island of stability in and around the 15-day orbit. However, direct integration of test 
particles showed that even orbits that were initially very stable in that region became 
unstable in later times. 

We would like to emphasize that although when comparing to Jovian-type planets, 
a test particle is a good approximation for a terrestrial-class object, the dynamical 
state of a test particle may not be a true representative of the dynamics of an actual 
Earth-sized body. The mutual interactions between this object and other bodies in 
the system play an important role in its dynamical state and may result in a stable 
orbital configuration in systems where test particle simulations indicate instability. In 
the case of GJ 876, however, we believe that the possibility of the existence of such an 
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0.50 1.00 1.50 2.00 

a[AU] 

Fig. 5 Graph of the stabihty of test particles in the region exterior to planet e using a grid 
of 1200 X 100. The color-coding represents the mean value of particles' MEGNO as described 
in Fig. [T] The nominal positions of the outer planets b and e are shown by white crosses. The 
white lines mark boundaries of crossing orbits with the three outer planets. Integrations were 
carried out for ~ 3400 years using a step-size of r = 0.1 days. The initial inclinations of test 
particles were set to i = 59° (i.e., co-planar with the orbits of the four planets). All the other 
orbital angles were set to zero. The initial orbital elements of the four planets were taken from 
Table H 

Earth-sized planet is very small. Given the current masses of the three giant planets 
of the system and their orbital configuration, in order for a fifth planet in the 15- 
day orbit to be stable, the mass of this object has to be larger than Earth-mass so 
that its mutual interactions with other planets would become significant to allow its 
(in)stability. Given the long history of the observation of GJ 876 (over 12 years), such 
a planet is expected to have been discovered by now. Results of our study suggest that 
the planetary system of GJ 876 is most likely dynamically full. 

It is also worth noting that from Hill radius calculations, one may argue that an 
additional, not too massive planet with a small eccentricity can maintain stability in 
the 15-day orbit. Although this argument has been shown to be valid in many multiple 
planet system, our numerical integrations indicate that in the system of GJ 876, Hill 
radius analysis may not be a good approach for identifying stability. The latter implies 
that this widely used technique has to be used with caution as it may not yield correct 
results at all times. 

The orbital architecture of the planets around GJ 876, with a super-Earth in a 
close-in orbit and three planets in a Laplace resonance, combined with the fact that this 
system is dynamically full, raises an important and interesting question; How did this 
planetary system form? Given that GJ 876 is an M star, it is very unlik ely that its outer 
three planets were formed in their current orbits. In fact, as shown bv lLauehlin et al~l 
( 20041 ) . the core-accretion model of giant planet formation fails to produce Jovian- type 
planets around M stars. The fact that many M stars host giant planets suggests that 
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these planets must have formed in ou ter regions of t he cir cums tellar disk and mig rated 
to their current orbits. As shown bv iLee fc Peale I ( 2002l l and lKlev et al. I (|2005l ). mi- 
grating planets may capture other bodies in their paths into mean-motion resonances 
(mainly the 1:2 resonance) where the two (or multiple) resonant bodies continue their 
migration until they either collide with the central star or reside in a stable orbit. For 
a more detailed and comprehensive review on the formation, mi gration, and resonanc e 
trapping of giant planets around M stars, we refer the reader to 'Haghighipour J ( 201ll l. 
While migrating inward, the resonant planets excite the orbits of smaller bodies caus- 
ing many of them to be scattered out of the system. Some of these bodies may also 
collide with one another and form larger objects (e.g., super-Earths). The interactions 
between these objects and the migrating planets may result in their scattering to stable 
close-in orbits. 
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